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Abstract Since its introduction by Gauss, Matrix Algebra has facilitated 
understanding of scientific problems, hiding distracting details and finding 
more elegant and efficient ways of computational solving. Todays largest prob- 
lems, which often originate from multidimensional data, might profit from even 
higher levels of abstraction. We developed a framework for solving tensor 
structured problems with tensor algebra that unifies concepts from ten- 
sor analysis, multilinear algebra and multidimensional signal processing. In 
contrast to the conventional matrix approach, it allows the formulation of 
multidimensional problems, in a multidimensional way, preserving structure 
and data coherence; and the implementation of automated optimizations of 
solving algorithms, based on the commutativity of all tensor operations. Its 
ability to handle large scientific tasks is showcased by a real- world, 4D medical 
imaging problem, with more than 30 million unknown parameters solved on a 
current, inexpensive hardware. This significantly surpassed the best published 
matrix-based approach. 

Keywords Tensor • Multidimensional problems • Tensor computations • 
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1 Introduction 

Computational problems with tensor structure, which involve large amounts 
of multidimensional data, arise in many fields of science and engineering [31[51 
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[2I1E1IIH1II1] ■ The standard way of dealing with such problems is not intrinsi- 
cally multidimensional, and assumes reduction of problem formulation to the 
classical matrix formalism, which is built on two-dimensional (matrix) and 
one-dimensional (vector) data structures. This is achieved by the reordering 
of multidimensional terms from the original problem formulation, into matrices 
and vectors. The derived matrix formulation is then treated using well-known 
techniques and algorithms of matrix algebra. Finally the solution is reconsti- 
tuted into its natural multidimensional form. 

While this approach reduces a multidimensional problem to the well-known 
standard form of matrix algebra, it introduces certain limitations. The result- 
ing formulation is no longer explicitly multidimensional, - a fact that may 
create difficulty in identifying and understanding important properties inher- 
ent to the particular problem [50]. Vectorization of multidimensional objects 
may lead to loss of spatial data coherence [SlfTB] , which can adversely affect the 
performance of solving algorithms. In many cases, the derived problem formu- 
lation does not have a straightforward, intuitive connection with the process 
of generating efficient solving algorithms. 

Currently, a growing interest in objects with more then two dimensions - 
tensors, has become evident in the scientific community. Introduced in Tensor 
Analysis and Multilinear Algebra, tensors gained the attention of practitioners 
of diverse fields (see the excellent review by Kolda [E]). However, translating 
tensor mathematics to a convenient computational framework raises many 
issues [191 ; including, convenient notation, ease of algorithm implementation, 
performance. 

Here we propose a comprehensive framework for solving tensor struc- 
tured problems by tensor algebra that allows natural and elegant for- 
mulation of multidimensional problems using multidimensional data structure 
directly. The proposed framework is based on a generalization of the concepts 
of matrix algebra to multiple dimensions; it incorporates and unifies existing 
approaches from multidimensional signal processing |5llT0l[2| , recent develop- 
ments in the field of tensor analysis [12] , and multilinear algebra [STlfTTlfTBllTTl 
[T^ITlfT^. It is based on ordered vector spaces and a generalized definition of 
tensor multiplication, based on the concept of tensor equations and the ten- 
sor inverse, and induces the design of fast computational solving algorithms 
with built-in expression analysis for automated generation of efficient imple- 
mentations on serial and parallel computers. In our companion paper we show 
how to computationally solve large real-world problems using our framework, 
and claim that this novel approach offers a natural, higher level abstraction 
to solving a broad range of very large computational problems which arise in 
multidimensional signal processing and other scientific fields. 



2 Tensor nomenclature, basic tensor operations 

In our framework, multidimensional data and multidimensional transforma- 
tions are described based on a formalism originating in Tensor Analysis [12] 
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and Physics [7] . Each object in the formulation represents a multidimensional 
entity - a tensor. Tensors describing multidimensional data belong to a space 
that is the outer product of a number of vector spaces fsee lA.3() . A tensor can 
be represented by a multidimensional array of its components^. For example, a 
scalar field defined on a discrete uniform three-dimensional grid with extents 
[NxtNy,Nz] can be expressed as a third order tensor 7"^^^, where tensor 
order is determined by the number of indices. Note that in our tensor nota- 
tion, direct correspondence between designation of tensor indices (x, y, z) and 
dimensions of the physical space {X, Y, Z) is typically used. Tensors indices 
can be either covariant (subscripts) or contravariant (superscripts), depending 
on the way the tensor components transform in respect to a change of basis 
(see lA.l)) . This property of a tensor index is called variance. For the common 
case of Euclidean spaces with orthonormal bases, there is no difference be- 
tween covariant and contravariant indices; thus the choice of index variance is 
influenced by the context and consistency of tensor expressions. 

Consider a transformation applied to tensor T^y^ which results in a tensor 
in the same tensor product space. The transformation can be expressed by 
•^xyz"^^ ■ Its application to T implies the use of classical tensor operations, 
such as the outer product, which expands tensor order 



AX-iy-iz-i ^xyz _ nx-iyizixyz (0 ^\ 

•^xyz ' ~ '^xyz K^-'^J 

and the contraction, which reduces tensor order 

QXiyiZixyz _ ^ ^ ^ ' ^ ^ QXiyiZixyz _ jyxiyizi ^2 2) 

X y z 

For designation of the contraction, we use the Einstein convention j7j which 
assumes implicit summation over a pair of indices with equal designation but 
different variance. From the practical point of view, it is convenient to combine 
the outer product and the contraction into a single operation, in analogy with 
the matrix multiplication. Thus by allowing simultaneous contraction over 
multiple corresponding index pairs we get 



j^xiyizi ^ j-xyz _ jyx^yizi ^2 3) 

Our framework relies on the fact that a tensor is an object which can 
implicitly contain all its properties as indices and component values. Due to 
this object nature in a practical implementation (|2.3p can be reduced to D := 
A*T, where "*" denotes tensor multiplication. 



2.1 Concept of ordered spaces, commutativity of tensor multiplication 

In some cases, multidimensional transformations can be decomposed to a prod- 
uct of one-dimensional transformations. Such transformations are called sepa- 



Further in the document the term tensor is identified with tensor components 
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rable. They are reduced to sequential application of one-dimensional transfor- 
mations from the decomposition. Using tensor formalism in three dimensions 
this can be expressed as 

j^xiyizi ^ j-xyz _ j^xi ^ Qyi ^ jyzi ^ q-xyz _ gxiyizi ^2 4) 

In conventional tensor formalism (and likewise in matrix formalism), the 
result of expression (j2.4p depends on the order of multiplicands (lack of com- 
mutativity). This is due to concatenation of non contracted indices, which is 
used for assurance of a unique result of the tensor product ^. In contrast, 
separable data handling operations, like filtering or resampling in signal pro- 
cessing, lead to the same result independently of the order in which separable 
transformations are applied. To avoid this formal, but practically important 
mismatch, and to add more flexibility to handling of tensor expressions, we in- 
troduce the following constraint: in analogy to physical space, which is ordered 
(right hand rule) , we require vector spaces used for construction of the tensor 
product to have a unique predefined order. This constraint leads to a funda- 
mental new property of the framework, compared to conventional approaches 
for matrix and tensors. The tensor multiplication becomes commuta- 
tive^ in addition to being associative. Formally this means that indices of the 
resulted tensor have unique positions determined by a predefined order of vec- 
tor spaces. Thus complex tensor products can be evaluated in arbitrary order 
but lead to the same result. For example, for ordered spaces X, F, Z we can 
write 

J\^x-i j^yi QZi ^q-xyz _ gyi .J^xi qzi ^q-xyz _ gyi qzi_ .j^xi _q-xyz _ gxiy-iz-i ^2 5) 

The commutativity of the tensor multiplication, achieved hereby, consid- 
erably increases the ease of handling tensor expressions. One of the important 
practical values of this property consists in simplification of semantic analysis 
of tensor expressions performed for the purpose of reduction of computational 
and storage requirements. For example, expression A^^* -B^^ -C^^ -Tt^^, [Nx < 
Ny < Nz < Nt, NxNy > Nt), which according to conventional tensor and 
matrix formalisms can be evaluated only in a single chain order, in our frame- 
work can be computed by reordering C^^ ■ {B^^ ■ (-4^^* • Tt^'^))- This gives 
the best computational performance with minimal memory requirements for 
storing intermediate results. At this point, it is possible to use the tensor mul- 
tiplication, in combination with tensor addition and subtraction, defined from 
multilinearity of tensor product space (see IA.4I) for formulating multidimen- 
sional problems in a natural and elegant multidimensional representation. Note 

^ With operations on mixed tensors, it can be rather difficult to track the overall order of 
array indices in the results. Thus in some tensor literature a special notation is used, where 
unique index order is defined by use of "dot" symbol {E'^y). But this brings an undesirable 
side-effect in significantly decreasing the readability of tensor expressions. 

^ Note that this constraint does not limit the expressiveness of tensors, as the order of 
vector spaces can be changed if necessary 



Computational Tensor Algebra 



5 



that formulations expressed in tensor terms are also convenient for differentia- 
tion in respect to multidimensional terms (see lA.8[) . an important prerequisite 
for multidimensional optimization problems. 

2.2 Extension of the tensor notation 

According to Einstein notation it is not permitted to have more than one 
index with same designation and variance. However, in our settings we do 
allow this, with agreement on some consistency rules. As an example, consider 
the following expression 

• = (2.6) 

This operation is a tensor analogue of the Khatri-Rao product 'SE^I which 
is a matching elementwise ( over i ) Kronecker product of two sets of vectors. 
Note that two formally equivalent indices are merged into a single one, without 
dependency on outer product or contraction applied to other indices. But the 
constraint is, that indices can be contracted only once after all possible merges 
resulting in a single pair of superscript and subscript. Here is an example 
which represents the most general form of the tensor product, combining outer 
product, contraction and elementwise product 

• • Bi ■ V,. ■ = {A\ • s;) • {v't ■ V,.. ■ w,) = niy ■ Mt. = 7^,. (2.7) 

Notice how neatly the proposed notation unifies different types of products 
which exist separately in Matrix Algebra (Matrix product, Kronecker product, 
Khatri-Rao product). 

3 Tensor equations and tensor inverse 

Expression A^^^^^'^ ■ U'-^y^' = g^iwi^i represents a relation between a multidi- 
mensional input U^y^- and a multidimensional output ^^iwi^i. This relation is 
determined by a transformation A%]y}^'^ . When the input is unknown, the ex- 
pression describes an inverse problem. In cases where a unique solution exists, 
the existence of a tensor inverse is implied - a transformation which maps A 
to the identity tensor 

xiyizi xyz xyz x\y\z\ \ / 

where = ■ Sy^ ■ is the identity transformation in three dimen- 

sions expressed in terms of Kronecker delta (see lA.Sp . Note that equation \3.1\ 
can be represented in an equivalent form as a matrix mapping from a vector 
space to a vector space, with proper reshaping of tensor terms. In cases when 
system tensor A has some structure, treatment of the problem using the tensor 
representation can be more advantageous than using the conventional matrix 
formalism, that is explained in the next sections. 
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4 Tensor solvers 

4.1 Direct tensor solvers, iterative tensor solvers, Krylov subspace tensor 
solvers 

Here we introduce tensor solvers which can be explicitly applied to solve such 
tensor equations computationally. The most straightforward way to find the 
solution tensor and, if necessary, a tensor inverse - is based on use of a tensor 
extension of Gauss elimination or a tensor analogue of LU decomposition - 
both belonging to the class of direct tensor solvers. These methods can benefit 
from the explicit tensor structure of equation coefficients, which is typically 
manifested by handling multidimensional sub-blocks of nonzero coefficients. In 
many practical problems, however, tensor equations are so large that solving 
them using direct methods becomes prohibitive. In this instance, we propose 
the class of iterative tensor solvers. For example, one can use a tensor Jacobi 
solver which is based on the iterative computation of tensor multiplication 
and belongs to the class of stationary iterative solvers. The pseudocode of the 
algorithm and an actual software implementation are shown on Fig. 14.11 

Note that the presented object based implementation of the Jacobi solver 
does not use explicit indexing of tensor components. All information about 
indices is contained in tensor objects which are defined once at the construction 
stage. This is sufficient for performing further computations with automatic 
semantic analysis for optimizing the performance, and without complicating 
the solving process by explicit use of indices. 

From the point of view of computational and storage requirements, the 
presented iterative solver is particularly attractive for application to tensor 
structured problems. As a basic example, consider the Poisson equation in 
three dimensions V'^u{x,y, z) — f{x,y,z) which after discretization on a uni- 
form grid can be expressed by the following tensor equation 

Here U is an unknown tensor, J- is the observation tensor, B, C are 
one-dimensional finite difference transformations which are equivalent to fil- 
ters with impulse response dependent on the order of the Laplacian approx- 
imation. This decomposition allows very efficient computation of the tensor 
multiplication by performing sequential separable transformations along indi- 
vidual dimensions of the tensor U. The separability renders the computation of 
the whole set of equation coefficients unnecessary, thus dramatically reducing 
storage requirements of the presented iterative algorithm (obviously, for such 
simple example, exploiting the tensor characteristics of the problem is also 
possible using the standard matrix approach with careful index bookkeeping) . 

Another type of iterative algorithms which are frequently used for solv- 
ing large inverse problems, is the family of Krylov subspace solvers, which 
includes GMRES, CG, BICGSTAB and other solvers [23]. In our extension 
of this solver class to tensors, most of them can be expressed as an iterative 
application of the tensor multiplication (see Fig. 14.21) and thus profit from the 
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for 1=0 until convergence 

y/xyz _ j/xyz _ qjx-iyizi . cxyz 

iA —LA ;V '^xiyizi 

'J^xiyizi _ jx\y\zx ^xyz _ j^xiyizi 

end 



(a) 



PROCEDURE TensorJacobi(A,B: Tensor; 

threshold: Scalar 
): Tensor; 

VAR R, E, U: Tensor; 
BEGIN 

E := InvMainDiag(A,B); 
R := A*U - B; 

WHILE R+*R > threshold DO 

U := U - R*E; 
R := A*U - B; 
END; 

RETURN U; 
END TensorJacobi; 



(b) 



(* define a world with vector spa<;es *) 



VAR 

w: World; 

A, B, C: Tensor; 
BEGIN 

NEW(w); 

w.DefineSpace('X',128) 
w.DefineSpace('Y',129) 
w.DefineSpace('Z',130) 

(* define system tensor *) 
A := Tensors.Laplacian(w,"x''l,x_,y'l,y_,z'l,z_"); 

(* define right hand side *) 
NEW(B,w,"x"l,y"l,z"l"); 

B. FloadC'rhs.dat"); 

(* solve the problem *) 
C := TensorJacobi(A,B,1.0E-4); 

END 



(c) 



Fig. 4.1 (a) Pseudocode of tensor Jacobi iterator for three dimensions {£ is the inverse 
of diagonal tensor constructed from main diagonal of .4 ). (b) Actual implementation of 

solver for any number of dimensions using our tensor library designed in Tensor Oberon 
programming language, which offers algebraic operators and optimized implementations for 
the basic tensor operations. Here "*" and "+*" stand for tensor multiplication and inner 
product respectively, (c) A simple code example of using the solver. 
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same benefits discussed above for the tensor Jacobi iterator. In contrast to 
Jacobi solver which only converges well for a limited class of problems, Krylov 
solvers perform significantly better for the large majority of problems. 

Fig. 14.21 presents a pseudocode and Oberon implementation of the tensor 
Conjugate Gradients (CG) algorithm applied to solve a 3D problem A%\Ji}^^ ■ 

Qxyz _ ^xiyizi ^ 



•pa^lMlzi 



^^h!'^' ■C^2'^ % initial residual 



: T^ajiaizi o/^ initial search directions 



for i=l by 1 until convergence 

Qxiyizi _ A^lVl^l . Cpxiyizi . I^ij/z ^ 
=ii — ^xyz V' ^12/121/ 

a 



% update solution approximation 

(^xyz _ Qxyz _j_ q, . ■pxiyizi . S^^y^zi 

Tlxiyizi ^ fixiyizi _ ^ . Qxiy\z\ % update residual 
% update search directions 

■pxiyizi _ -J^xiyizi _|_ / [n \ . -pxiyizi 



p ■■ 

end 



Pl 



R 
P 
D 

rho 



: B-A*C; 
R; 

: Delta([C.inds, -B.inds]); 
= R+*R; 



REPEAT 

Q := A*(P*D); 
alpha := rho/(P+*Q); 

C := C + alpha*P*D; 
R := R - alpha*Q; 
rhol = R+*R; 

P := R + (rhol/rho)*P; 

rho := rhol; 
UNTIL rho > threshold 



Fig. 4.2 Pseudocode and given side by side actual implementation of Conjugate Gradient 
iterator 



4.2 Tensor Multigrid algorithms unify algebraic multigrid and geometric 
multiresolution 

In the matrix domain, multigrid methods for solving linear system of equa- 
tions are known for their fast convergence and computational efficiency. These 
methods perform iterative improvement of a solution approximation, based 
on smoothing of the solution error on multiple scales of the problem. Basi- 
cally, there are two types of multigrid/multiresolution strategies: geometric 
multiresolution preserves spatial coherence in multidimensions by working on 
datasets at different resolutions but its formulation is often problem specific. 
Algebraic matrix multigrid algorithms (AMG) offer a 'black box' approach 
for solving matrix equations, but they rely on values of matrix coefficients 
only and risk neglecting data coherence, at a potential cost of degraded per- 
formance. 

We have used tensor equations to express both approaches in a unified 
fashion. The proposed tensor multigrid algorithm (TMG) has an important 
advantage over AMG: in contrast to the matrix algorithm, which neglects 
multidimensionality and the structure of the problem, TMG unifies both con- 
cepts of geometric and algebraic multiresolution, by preserving and exploiting 
multidimensionality and spatial data coherence (see Fig. 14. 3p . 
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Fig. 14.31 shows pseudocode for a Tensor Algebraic Multigrid V-cycle algo- 
rithm applied to solve a 3D tensor problem -4i^y|^^^ .Q'^v^ = i^xivi^i^ Problems 
in higher dimensions are handled in a similar way. Fig. 14.41 compares visually 
how AMG and TMG are applied to a problem of reconstruction of a two- 
dimensional image from a few arbitrarily taken samples, using non- uniform 
spline interpolation [3 . This example shows that performance of the AMG 
algorithm is significantly degraded due to data distortion, which is introduced 
at coarse matrix scales, whilst the TMG converges to a good approximation of 
the solution within a few iterations, due to its preservation of spatial coherence. 



function C^^ = TMGVCycle(scaZe,yl^i^i^i ,C^!'^,B^i!'i^i) 
if scale ^ last then 

(^xyz ^ Presmooth(^;;:lIl^^C^!'^B^l^'l^l) 

% get preliminary constructed system scale and 
% corresponding restriction/prolongation operators 

[a112^^ , C"--, B"!"!™! , V^l , K,S^l , Sy,Tr,' , f£\ = GetScale(scaie) 

% separable restriction of the residual 

^uvw ^ TMGVCycle(sca/e-|- l.^^i^i^SC"'"",^"!''!™!) % V-cycle on next scale 
^xyz ^ (.xyz ^x , gy , fz , ^uvw % ^^^^^^ ^^-^ correction 

Qxyz ^ PostsmOOth(^^lIl^^C^!'^B^l!'l^l) 

else 

C^f^ = DirectSolveC^Jili^SB^ifi^i) 
end 
end 

Fig. 4.3 Pseudocode of a TMG V-cycle algorithm for three dimensions 



5 Automatic tensor expression analysis for high-performance 
implementations 

Because the high level tensor formulation contains the complete tensorial struc- 
ture information of a complex mathematical problem, automatic expression 
analysis can be performed, resulting in automated software optimization and 
code generation. We have experimentally verified this aspect in our tensor 
library and have found that such automatic analysis is capable of optimiz- 
ing code, with regard to computational and memory requirements, in both 
a problem-specific and hardware-specific manner. This result profits strongly 
from the commutativity of all tensor operations the significant differentiating 
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Projection from solution 2 iterations 15 iterations iterations 

at the coarsest scale 



Fig. 4.4 Performance of matrix (AMG) and tensor (TMG) multigrid algorithms in appli- 
cation to a problem of image reconstruction from an incomplete set of points. The tensor 
approach (TMG), through preservation of data locality in the coarser scales, yields higher 
quality solutions and faster convergence. 

feature over standard matrix formulations which are non-commutative with 
respect to matrix multiplication. Some aspects of such automated code gener- 
ation have also been discussed in . This particular topic is a promising area 
for future research. 



6 Application of the framework to large real-world problems 

Our framework was successfully applied to solve a large real-world computa- 
tional problem - a spline-based global variational reconstruction [3] of a mul- 
tidimensional signal from incomplete and spatially scattered measurements in 
four dimensions. The problem was originally formulated in terms of matrix al- 
gebra. In four dimensions, the matrices involved in the computations become 
extremely large, essentially prohibiting computation of the problem on current 
standard hardware. For example, for a moderate size of the reconstruction grid 
128 X 128 X 128 x 16 we have to deal with a matrix which is represented by 
7'514'144'528 nonzero elements (in compressed block-band diagonal storage, 
30 GBytes in single precision). Using our framework, we reformulated the prob- 
lem in terms of tensors. The resulting tensor formulation helped to analyze 
the mathematical structure of the problem and to derive its decomposition 
in the form of a 1-rank tensor decomposition known as Canonical Decom- 
position (CANDECOMP) [TTlfMllTlfTQ] . More details about proposed tensor 
formulation and solving algorithm can be found in our companion paper. 

The identified tensorial structure allowed efficient computation of the ten- 
sor multiplication which therefore can be efficiently used within a tensor Krylov 
solver. Another important implication of the identified property of the tensor 
formulation, is a dramatic reduction of storage requirements. Our approach 
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does not require explicit storage either of the system tensor or its decompo- 
sition, and thus the problem of size, discussed above (33'554'432 unknowns 
parameters), for about 9'000'000 measurements, can be computed on current 
inexpensive multi core computer equipped with only 2 GBytes of physical 
memory in less than 60 minutes, surpassing the capability of a published, 
matrix-based solving algorithm [3 for the same problem by far. The decom- 
posability of the tensor multiplication allows convenient and complete paral- 
lelization of computations using all currently available parallelism paradigms, 
including Single Instruction Multiple Data (SIMD), Multi Core, Clusters of 
PCs and General Purpose GPU (GPGPU) computing technology. 

Note that parallelizability of tensor operations, in general, has proven to 
be very beneficial: on Multi Core CPUs we have observed a speed increase 
almost proportional to the number of cores; and on a recent graphics card (ATI 
Radeon HD 4870 X2) an additional increase has been measured for specific 
examples, of more than 10, compared to a Quad Core CPU. 



7 Discussion 



In this work we introduced a framework for computational solving of large ten- 
sor structured problems. The proposed approach leads to a natural, dimensionality- 
preserving formulation of tensor structured problems, and by maintaining the 
multidimensional structure of the data, and the commutativity properties of 
the framework, directly induces computationally efficient solving strategies 
that can profit from automatic expression analysis, separability properties of 
the tensor formulation, and the preservation of spatial coherence of the data, 
to speed up convergence in tensor solvers as compared to their matrix coun- 
terparts. 

The chosen formalism based on Einstein notation together with introduced 
commutativity of tensor multiplication, makes the framework well suited for 
implementing a tensor programming language, offering self-optimized compu- 
tations of tensor expressions by semantic analysis of their terms while avoiding 
explicit use of tensor indexing in the actual implementations of solving pro- 
grams. This was verified by implementing our own tensor library built on the 
principles of the proposed framework. 

The properties presented above distinguish our framework from one often 
used in the field of tensor decompositions (271[T71[TB1[TT1[T^ [T1[T^ . but do not 
limit its use for solving problems studied in that field. 

The ongoing research shows that the framework is well suited for solving 
problems of signal processing in higher dimensions, which have inherent tensor 
structure. We are currently investigating its usefulness in problems from image 
reconstruction, computer vision and fluid dynamics. 
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In summary, this high-level approach fits well into the landscape of increas- 
ing interest in tensors as a workhorse for elegant and efficient solving of those 
scientific problems which arise from multidimensional data. 



Table 7.1 Basic multidimensional data processing building blocks in separable tensor form 



Signal processing operator 


Tensor expression 


Corresponding 
implementation ^ 


Discrete convolution 


£,Q,'H - ID convolution transforms 


F := E*G*H*C; 


Finite difference 


jrxyiz _ gyi . ^xyz 

jrxyz\ ^hi^^ ■ (^-^y^ 

CH - ID finite difference transforms 


F := E*C; 
F := G*C; 
F := H*C; 


Discrete Fourier transform 


j^^ivi^i = ^•^i . gyi . -^^1 . Q'-^yz 
£,g,H- ID DFT transforms 


F := E*G*H*C; 


Rotation 


jrx2yizi ^ 2?J2 . ^^^1 . gyi . -^^i) . (jaii/z 
f , £,Q,'H - ID shear transforms 


F := D*E*G*H*C; 


Upsampling/ downsampling 


■puvw _ gv _ . Qxyz 

£,0,^-1 - ID up-/downsampling transforms 


F := E*G*H*C; 



Table 7.2 Large scientific problems suited for solution by computational Tensor Algebra 



Computational Application 


Scientific Field 


Multilinear data models using tensor decompositions 
(e.g. CANDECOMP, TUCKER, HOSVD): 

- Pattern recoenition 1240281 

- Data compression 1131 

- Independent Component Analysis [6lll5| 

- Signal filtering [20] 

- Modeling of fiuorescence data ^ 

- Social Network Analysis [l] 


Signal Processing 
Machine Learning 
Computer Vision 
Data Mining 
Chemistry 
Neuroscience 


Solving inverse multidimensional problems: 

- Least squares problems 3*21] 

- Differential equations (Poisson, Navier-Stokes, Finite 
Element Methods) 

- Integral equations [8] 


Signal Processing 
Computer Vision 
Machine Learning 
Fluid Dynamics 
Chemistry 
Electromagnetism 


Modeling of properties of complex systems consisting of many 
interacting elements expressed by tensor expressions: 
- Modeling of electronic and optical properties of 
molecules and their interactions |25II22| 


Computational Physics 
Quantum Chemistry 
Computational Chemistry 
Material Science 



* Objects C,D,E,F,G,H correspond respectively to C, "D, f , J^, (7, W from tensor expres- 
sions; "*" denotes the tensor product 

^ The tensor representation of finite difference is the base for derivative related operators 
such as gradient, divergence, curl and Laplacian. 
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A APPENDIX 

A.l Dual vectors spaces, contravariant and covariant mechanisms of 
transformation 

Let y be a real vector space with finite number of dimensions Nv ■ We designate by V* 
unique vector space of same dimensionality that is dual in respect to V . Any element t* G V* 
is a linear map V — )■ 5R 

(t*,t> =t*(t) e 5R (A.l) 

For a given base of V transformations to a new base e^-^ and vice versa are defined 

by 

Nv Nv 

where and S^^ ^.re direct and indirect base transformations respectively. A vector 

from V given by Y.v=i ■ according to llX2t in a new coordinate system is represented 
by components 

Nv 

r"i = ^ ■ (A.s) 

u = l 

Components of such vectors which transform indirectly in respect to transformation of 
the base are called contravariant. 

For the base of V and dual base e"i of V* the following equation is satisfied 

(e-,e„).c ={;;;;:;;: (a.4) 

where <5JJ^ is Kronecker delta. Thus the dual base e"i transforms according to 
Nv Nv 

vi=\ 112 = 1 

In contrast to contravariant vector components, components of vectors from dual space 
V* transform in the same way as base of V 

Nv 

r., = Yl ^2 ■ Til (A-6) 

VI —1 

Components of such vectors are called covariant. 



A.2 Generalized definition of a tensor 

Definition A.l Let ordered set of real vector spaces Ui, . . . ,Up , Vi, . . . , Vq with finite 
dimensions I\, . . . ,Ip , Ji, . . . , Jq and their dual vector space l/J', . . . , Up , , . . . , Vq be 
given. By identifying (C/*)* with Ui we define every vector g Ui to be a linear map 
U* — > 5i given by 

u,(u*) = (u*,u,) G R (A.7) 

for arbitrary u* g U* , where (, ) denotes the inner product. With (P + Q) vectors 
ui e l/i , . . . , Up G Up ,v1&V*,..., Vq G Vq let the element denoted T = ui X ■ ■ ■ X up X 
Vj X ■ ■ ■ X Vq be a (P + (3)-linear map from l/i X ■ • • X Up X Vj* X ■ ■ ■ X Vq to SR defined by 
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ui X ■ ■ ■ X Up X V* X ■ ■ ■ X Vq (a*, . . . , aj,, bi, . . . , bg) = 
(aj,ui) ■ ■ ■ (ap,up)(vf,bi) ■ ■ ■ (v^.b^) 

where a* , hj are arbitrary vectors in U* and Vj respectively. The element T is termed a 
(P+Q)-th order decomposed tensor, P-times contravariant and Q-times covariant. The space 
generated by all linear combinations of decomposed tensors is termed the tensor product 
space. Any arbitrary tensor can be represented as a weighted sum of decomposed tensors 
[TTl[Tlll9yi5l . This definition is compatible with, and extends [T7l[2T) . 



A. 3 Tensor notation 



To each dimension of the physical space we assign a vector space with a given dimensionality. 
For example for 3-dimensional physical space we define vector spaces X, Y, Z with sizes 
-^'^x 1 1 ■ In such vector spaces one can define multiple coordinate systems related to 
each other by transformations discussed above. Components of a tensor in some tensor 
product space are designated as an array with multiple indices. Each index in the designation 
has direct correspondence to a dimension of the physical space and thus to one of the defined 
vector spaces. For example for 3-dimensional case a tensor which is an element oi X xY X Z 
can be designated by ^^iVj'^k Subindices i, j, k are used for distinguishing between different 
coordinate systems. 



A. 4 Elementwise tensor operations 

Linearity of the defined tensor product space implies definition of the following elementwise 
tensor operations such as tensor addition, subtraction and multiplication by a scalar; 



■^xyz ' ^xyz ^xyz \^^"-^J 
•^xyz '~^xyz ^xyz 

\ . A^lVl^l — (-xiyizi CA iTl 

^ -^xyz *^xyz \^-^^) 

With introduced constraint of ordered vector spaces all properties of these operations 
such as associativity and commutativity remain unchanged. 



A. 5 Kronecker delta 

Kronecker delta presented in expression llA.411 is a second order tensor which is equivalent 
to the identity transformation. Mixed version of the tensor applied to components of a 
contravariant or covariant first order tensor (vector) does not introduce changes 



5%^ ■ = r^i = (A.12) 

■Ux,=Ux=Ux^ (A.13) 

Covariant or contravariant versions of the Kronecker delta change variance of tensor 
indices, but do not change values of tensor components 



■Ux^ 



— ^X 



(A.14) 
(A.15) 
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Note that here indices x and xi correspond to identical coordinate systems. In the same 
way Kronecker delta can be applied along some dimension of a higher order tensor 

■ T'"'^^ = ■ T'^y^ ■ ■ T''^^ = T''^" (A.16) 
The tensor multiplication of Kronecker deltas for different vector spa<;es forms multidi- 
mensional identity transformation 



^xyz 
xyz 



. q~xyz 



(A.17) 
(A.18) 



A. 6 Inner product 

the inner product of two third order Euclidean tensors and B^^^ is defined by 

{A, B) = {A-y-^ ■ 5,,^,yy,,,,) ■ (B-^- . -1) = i-is'i-i . B.,y,., (A.19) 

At the same time using the language of elementwise products the inner product caji be 
expressed as 

{A, B) = {A-y^ ■ B-y^) ■ {S^x^yy^zz^ ■ = C^l^l^l ■ X.i„zi (A.20) 

where Tx^y^z^ is a tensor with all components equal to one, which can be easily checked 
based on the properties of Kronecker delta. 

Note that in general case of non-Euclidean spaces the inner product is defined with the 
use of the metric tensor which for Euclidean case reduces to the Kronecker delta. 



A. 7 Tensor transposition 

In classical settings transposition of a tensor is defined by changing positions of tensor 
indices. But in our framework the order of indices is unique a<;cording to a predefined order 
of vector spaces. Thus patterns like A^A from matrix normal equations correspond to the 
tensor multiplication with change of variance by the Kronecker delta in Euclidean spa<;es or 
the metric tensor in general case. 



A. 8 Differentiation in respect to a tensor 

When dealing with an optimization problem, differentiation in respect to the optimization 
parameter may be required. Consider an example of diff'erentiating a tensor valued function 
of the form /(W) = ^^J^^^i • 



d_ 

where indices {x,y,z) and {x2,y2,Z2) correspond to identical coordinate systems. Using 
tensor transformational properties it is straightforward to show that 



dW2y2Z2 ^2y2Z2 

Thus we have 



gi^xyz 

&III2Z2 (A-22) 



_f(1l\ — A^IVIZI , xxyz _ AXiyizi _ jtxiyizi /a r,o\ 

J J — -^xyz X2y2Z2 '^X2y2Z2 — "^xyz 



d_ 

'au' 

In a similar way one can compute derivative of more complex scalar and tensor valued 
functions. 
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